/* This program is free software: you can redistribute it and/or modify
   it under the terms of the GNU General Public License as published by
   the Free Software Foundation, either version 3 of the License, or
   (at your option) any later version.

   This program is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.

   You should have received a copy of the GNU General Public License
   along with this program.  If not, see <http://www.gnu.org/licenses/>. */
#include <windows.h>
#include <process.h>
#include <math.h>
#include "Mandelbrot.h"
#include "Render.h"

/* MandelbrotPixel() performs a test on the pixel represented by the complex
   number c, iterating at most MAX_ITERATIONS times, and testing to see after 
   each iteration if the number c is diverging to infinity (i.e. isn't in the 
   set).  In reality, you just check to see if the modulus is bigger than some
   fixed value. The formula applied on each iterations is just: 
   z[n+1] = z[n]^2 c, z[0] = 0 */
COLORREF MandelbrotPixel(const COMPLEX c)
{
	UINT i, colour;
	COMPLEX prev;
	COMPLEX curr;
	const double scaleFactor = 255.0 / MAX_ITERATIONS;

	prev.imag = prev.real = 0.0;
	for (i = 0; i < MAX_ITERATIONS; i++)
	{
		/* z^2 - for the number (a + bi) this results in a^2 - b^2 + 2abi */
		curr.real = prev.real * prev.real - prev.imag * prev.imag;
		curr.imag = 2.0 * prev.real * prev.imag;
		/* + c */
		curr.real += c.real;
		curr.imag += c.imag;

		/* the bailout condition is if the modulus is larger than the
		   bailout value, we assume the sequence in diverging */
		if (ComplexModulus(curr) > BAILOUT_VALUE)
		{
			/* scale the colour (Red, Green or Blue value) uniformly across 
			   the range [0, MAX_ITERATIONS] */
			colour = (int) (255.0 - scaleFactor * i);
			return RGB(colour, colour, 255);
		}
		prev = curr;
	}
	return RGB(0, 0, 0);
}

/* NormaliseVal() takes a screen x or y coordinate, and normalises it to 
   the given interval, say [-1,2] or [-1,3].  The current screen resolution
   is also required, so that you know how much to scale by */
double NormaliseVal(const double dIntervalStart, const double dIntervalEnd,
                    const long lRes, const long lPos)
{
	double dIntervalWidth, dNormalised, dRes, dPos; 
	dRes = (double) lRes;
	dPos = (double) lPos;

	dIntervalWidth = dIntervalEnd - dIntervalStart;
	dNormalised = dPos / (dRes / dIntervalWidth);
	dNormalised += dIntervalStart;

	return dNormalised;
}

/* ComplexModulus() gets the modulus of a given complex number, i.e.
   sqrt(real^2 + imag^2).  However, the actual computation done is
   slightly differently to avoid overflow -
   |a + ib| = |a|sqrt(1 + (b/a)^2 if |a| >= |b|
   or         |b|sqrt(1 + (a/b)^2 if |b| < |a| */
double ComplexModulus(const COMPLEX complexVal)
{
	double modulus;
	const double a = complexVal.real, b = complexVal.imag;

	// avoid division by 0
	if (a == 0.0 && b == 0.0)
		return 0.0;
	else if (fabs(a) >= fabs(b))
		modulus = fabs(a) * sqrt(1 + pow((b/a), 2));
	else
		modulus = fabs(b) * sqrt(1 + pow((a/b), 2));
	
	return modulus;
}

/* MandelbrotRender() takes an area of the screen, and does all the 
   calculations required for this area.  Each pixel is converted 
   into a complex number and tested to see if it's in the set, with the colour
   representing how much the pixel has diverged from 0; black - in the set,
   white - not */
void MandelbrotRender(THREADPARAM *pParams)
{
	int xScreen, yScreen, xBmp, yBmp;
	double imagMax, imagMin, realMax, realMin;
	COMPLEX currPixel;

	imagMax = pParams->range.imagMax;
	imagMin = pParams->range.imagMin;
	realMax = pParams->range.realMax;
	realMin = pParams->range.realMin;

	pParams->threadStatus = RENDERING;
	for (xScreen = pParams->plotArea.left, xBmp = 0; 
	     xScreen < pParams->plotArea.right; xScreen++, xBmp++)
	{
		/* check to make sure the thread hasn't been cancelled */
			if (pParams->threadStatus == CANCELLED)
				break;

		/* convert from screen coordinate to cartesian plane, first finding the 
		   real component */
		currPixel.real = 
			NormaliseVal(realMin, realMax, pParams->xImageRes, xScreen);
		for (yScreen = pParams->plotArea.bottom, yBmp = 0; 
		     yBmp < pParams->plotArea.bottom; yScreen--, yBmp++)
		{
			/* find the imaginary component, the image is flipped vertically 
			   because the Windows coordinate system has 0 as the top */
			currPixel.imag = NormaliseVal(imagMin, imagMax, pParams->yImageRes,
			                              yScreen);
			
			SetPixel(pParams->hdcRender, xBmp, yBmp, 
			         MandelbrotPixel(currPixel));
		}
	}
}